#ifndef FUNCTION_HPP
#define FUNCTION_HPP

#include <vector>
#include <cmath>
#include <stdexcept>

/**
 * @brief Abstract base class for mathematical functions
 */
class Function {
public:
    virtual ~Function() = default;

    /**
     * @brief Pure virtual function to evaluate the function at a given point x
     * @param x Input point
     * @return Value of the function at x
     */
    virtual double operator()(double x) const = 0;

    /**
     * @brief Virtual function to evaluate the derivative of the function at a given point x
     * @param x Input point
     * @param h Step size for the finite difference method (default is 1e-7)
     * @param n Order of the derivative (default is 1)
     * @return Value of the derivative at x
     */
    virtual double derivative(double x, double h = 1e-7, int n = 1) const {
        if (n == 0) {
            return (*this)(x);
        } else if (n == 1) {
            return ( (*this)(x + h) - (*this)(x - h) ) / (2 * h); 
        } else {
            return (derivative(x + h, h, n - 1) - derivative(x - h, h, n - 1)) / (2 * h);
        }
    }
};

#endif // FUNCTION_HPP